/* CS580 Homework 4 */

#include	"stdafx.h"
#include	"stdio.h"
#include	"math.h"
#include	"Gz.h"
#include	"rend.h"


int GzRotXMat(float degree, GzMatrix mat)
{
// Create rotate matrix : rotate along x axis
// Pass back the matrix using mat value
	mat[1][1]=cos(degree*PI/180);
	mat[1][2]=-sin(degree*PI/180);
	mat[2][1]=sin(degree*PI/180);
	mat[2][2]=cos(degree*PI/180);

	return GZ_SUCCESS;
}


int GzRotYMat(float degree, GzMatrix mat)
{
// Create rotate matrix : rotate along y axis
// Pass back the matrix using mat value
	mat[0][0]=cos(degree*PI/180);
	mat[2][0]=-sin(degree*PI/180);
	mat[0][2]=sin(degree*PI/180);
	mat[2][2]=cos(degree*PI/180);

	return GZ_SUCCESS;
}


int GzRotZMat(float degree, GzMatrix mat)
{
// Create rotate matrix : rotate along z axis
// Pass back the matrix using mat value
	mat[0][0]=cos(degree*PI/180);
	mat[0][1]=-sin(degree*PI/180);
	mat[1][0]=sin(degree*PI/180);
	mat[1][1]=cos(degree*PI/180);

	return GZ_SUCCESS;
}


int GzTrxMat(GzCoord translate, GzMatrix mat)
{
// Create translation matrix
// Pass back the matrix using mat value
	mat[0][3]=translate[0];
	mat[1][3]=translate[1];
	mat[2][3]=translate[2];

	return GZ_SUCCESS;
}


int GzScaleMat(GzCoord scale, GzMatrix mat)
{
// Create scaling matrix
// Pass back the matrix using mat value
	mat[0][0]=scale[0];
	mat[1][1]=scale[1];
	mat[2][2]=scale[2];

	return GZ_SUCCESS;
}


//----------------------------------------------------------
// Begin main functions

int GzNewRender(GzRender **render, GzRenderClass renderClass, GzDisplay	*display)//
{
/*  
- malloc a renderer struct /
- keep closed until all inits are done /
- setup Xsp and anything only done once /
- span interpolator needs pointer to display /
- check for legal class GZ_Z_BUFFER_RENDER /
- init default camera /
*/ 
	int i,j,k;
	GzCoord x,y,z;
	float cl,norm;

	(*render) = (GzRender*) operator new (sizeof(GzRender));
	(*render)->renderClass = renderClass;
	(*render)->display = display;
	(*render)->interp_mode = 0;
	(*render)->matlevel = 0;
	(*render)->numlights = 0;
	(*render)->open = 0;				/*Keep closed*/
	(*render)->spec = DEFAULT_SPEC;
	(*render)->tex_fun = NULL;
	(*render)->aashiftx = 0;
	(*render)->aashifty = 0;

	(*render)->Ka[RED] = 0.1;			/*Throw DEFAULT_AMBIENT, DEFAULT_DIFFUSE, DEFAULT_SPECULAR to Ka, Kd, Ks*/
	(*render)->Ka[GREEN] = 0.1;
	(*render)->Ka[BLUE] = 0.1;
	(*render)->Kd[RED] = 0.7;
	(*render)->Kd[GREEN] = 0.6;
	(*render)->Kd[BLUE] = 0.5;
	(*render)->Ks[RED] = 0.2;
	(*render)->Ks[GREEN] = 0.3;
	(*render)->Ks[BLUE] = 0.4;

	/*Initialize camera by defualt*/
	(*render)->camera.FOV = DEFAULT_FOV;
	(*render)->camera.lookat[X] = 0.0;				/*Init lookat position*/
	(*render)->camera.lookat[Y] = 0.0;
	(*render)->camera.lookat[Z] = 0.0;
	(*render)->camera.position[X] = DEFAULT_IM_X;	/*Init world position*/
	(*render)->camera.position[Y] = DEFAULT_IM_Y;
	(*render)->camera.position[Z] = DEFAULT_IM_Z;
	(*render)->camera.worldup[X] = 0.0;				/*Init world up vector*/
	(*render)->camera.worldup[Y] = 1.0;
	(*render)->camera.worldup[Z] = 0.0;
	for(i=0;i<4;i++)								/*Init Xpi matrix*/
	{
		for(j=0;j<4;j++)
		{
			if(i==j)
			{
				(*render)->camera.Xpi[i][j] = 1.0;
			}
			else
			{
				(*render)->camera.Xpi[i][j] = 0.0;
			}
		}
	}
	(*render)->camera.Xpi[3][2] = tan(((*render)->camera.FOV/2)*PI/180);

	for(i=0;i<4;i++)								/*Init Xiw matrix*/
	{
		for(j=0;j<4;j++)
		{
			if(i==j)
			{
				(*render)->camera.Xiw[i][j] = 1.0;
			}
			else
			{
				(*render)->camera.Xiw[i][j] = 0.0;
			}
		}
	}

	/*According to prof, which is SideVec = ViewDir x WorldUp, and UpVector = SideVec x ViewDir*/
	cl = sqrt(((*render)->camera.lookat[X]-(*render)->camera.position[X])*((*render)->camera.lookat[X]-(*render)->camera.position[X])
					+((*render)->camera.lookat[Y]-(*render)->camera.position[Y])*((*render)->camera.lookat[Y]-(*render)->camera.position[Y])
					+((*render)->camera.lookat[Z]-(*render)->camera.position[Z])*((*render)->camera.lookat[Z]-(*render)->camera.position[Z]));
	z[X] = ((*render)->camera.lookat[X]-(*render)->camera.position[X])/cl;
	z[Y] = ((*render)->camera.lookat[Y]-(*render)->camera.position[Y])/cl;
	z[Z] = ((*render)->camera.lookat[Z]-(*render)->camera.position[Z])/cl;				/*Init Z*/
	norm = sqrt((z[Z]*(*render)->camera.worldup[Y]-z[Y]*(*render)->camera.worldup[Z])*(z[Z]*(*render)->camera.worldup[Y]-z[Y]*(*render)->camera.worldup[Z])
				+(z[X]*(*render)->camera.worldup[Z]-z[Z]*(*render)->camera.worldup[X])*(z[X]*(*render)->camera.worldup[Z]-z[Z]*(*render)->camera.worldup[X])
				+(z[Y]*(*render)->camera.worldup[X]-z[X]*(*render)->camera.worldup[Y])*(z[Y]*(*render)->camera.worldup[X]-z[X]*(*render)->camera.worldup[Y]));
	x[X] = -(z[Y]*(*render)->camera.worldup[Z]-z[Z]*(*render)->camera.worldup[Y])/norm;
	x[Y] = -(z[Z]*(*render)->camera.worldup[X]-z[X]*(*render)->camera.worldup[Z])/norm;
	x[Z] = -(z[X]*(*render)->camera.worldup[Y]-z[Y]*(*render)->camera.worldup[X])/norm;	/*Init X*/
	norm = sqrt((x[Z]*z[Y]-x[Y]*z[Z])*(x[Z]*z[Y]-x[Y]*z[Z])+(x[X]*z[Z]-x[Z]*z[X])*(x[X]*z[Z]-x[Z]*z[X])+(x[Y]*z[X]-x[X]*z[Y])*(x[Y]*z[X]-x[X]*z[Y]));
	y[X] = -(x[Y]*z[Z]-x[Z]*z[Y])/norm;
	y[Y] = -(x[Z]*z[X]-x[X]*z[Z])/norm;
	y[Z] = -(x[X]*z[Y]-x[Y]*z[X])/norm;													/*Init Y*/
	/*******************************************************************************************/

	/*Compute Xiw*/
	(*render)->camera.Xiw[0][0] = x[X];
	(*render)->camera.Xiw[0][1] = x[Y];
	(*render)->camera.Xiw[0][2] = x[Z];
	(*render)->camera.Xiw[1][0] = y[X];
	(*render)->camera.Xiw[1][1] = y[Y];
	(*render)->camera.Xiw[1][2] = y[Z];
	(*render)->camera.Xiw[2][0] = z[X];
	(*render)->camera.Xiw[2][1] = z[Y];
	(*render)->camera.Xiw[2][2] = z[Z];
	(*render)->camera.Xiw[0][3] = -(x[X]*(*render)->camera.position[X]+x[Y]*(*render)->camera.position[Y]+x[Z]*(*render)->camera.position[Z]);
	(*render)->camera.Xiw[1][3] = -(y[X]*(*render)->camera.position[X]+y[Y]*(*render)->camera.position[Y]+y[Z]*(*render)->camera.position[Z]);
	(*render)->camera.Xiw[2][3] = -(z[X]*(*render)->camera.position[X]+z[Y]*(*render)->camera.position[Y]+z[Z]*(*render)->camera.position[Z]);
	/******************************/

	/*Initialize color*/
	for(i=0;i<3;i++)
	{
		(*render)->flatcolor[i] = 0;
	}

	/*Initialize the lights*/
	for(i=0;i<MAX_LIGHTS;i++)
	{
		for(j=0;j<3;j++)
		{
			(*render)->lights[i].color[j] = 0.0;
			(*render)->lights[i].direction[j] = 0.0;
		}
	}

	/*Initialize the ambient light*/
	for(i=0;i<3;i++)
	{
		(*render)->ambientlight.direction[i] = 0.0;
		(*render)->ambientlight.color[i] = 0.0;
	}

	/*Initialize matrix stack*/
	for(i=0;i<MATLEVELS;i++)
	{
		for(j=0;j<4;j++)
		{
			for(k=0;k<4;k++)
			{
				if(k==j)
				{
					(*render)->Ximage[i][j][k] = 1.0;
					(*render)->Xnorm[i][j][k] = 1.0;
				}
				else
				{
					(*render)->Ximage[i][j][k] = 0.0;
					(*render)->Xnorm[i][j][k] = 0.0;
				}
			}
		}
	}

	/*Initialize Xsp matrix*/
	for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			(*render)->Xsp[i][j] = 0.0;
		}
	}
	(*render)->Xsp[0][0] = display->xres/2;
	(*render)->Xsp[0][3] = display->xres/2;
	(*render)->Xsp[1][1] = -display->yres/2;
	(*render)->Xsp[1][3] = display->yres/2;
	(*render)->Xsp[2][2] = INT_MAX*tan(((*render)->camera.FOV/2)*PI/180);
	(*render)->Xsp[3][3] = 1.0;
	/***********************/

	if(!GZ_Z_BUFFER_RENDER)
	{
		return GZ_FAILURE;
	}
	return GZ_SUCCESS;

}


int GzFreeRender(GzRender *render)
{
/* 
-free all renderer resources
*/
//	GzFreeDisplay(render->display);
	delete(render);
	return GZ_SUCCESS;
}


int GzBeginRender(GzRender *render)
{
/*  
- set up for start of each frame - clear frame buffer /
- compute Xiw and projection xform Xpi from camera definition /
- init Ximage - put Xsp at base of stack, push on Xpi and Xiw /
- now stack contains Xsw and app can push model Xforms if it want to. /
*/
	int i,j,k;
	GzCoord x,y,z;
	float cl,norm;

	GzInitDisplay(render->display);					/*Set up for start of each frame*/

	/*Compute Xiw and Xpi*/
	for(i=0;i<4;i++)								/*Init Xpi matrix*/
	{
		for(j=0;j<4;j++)
		{
			if(i==j)
			{
				render->camera.Xpi[i][j] = 1.0;
			}
			else
			{
				render->camera.Xpi[i][j] = 0.0;
			}
		}
	}
	render->camera.Xpi[3][2] = tan((render->camera.FOV/2)*PI/180);

	for(i=0;i<4;i++)								/*Init Xiw matrix*/
	{
		for(j=0;j<4;j++)
		{
			if(i==j)
			{
				render->camera.Xiw[i][j] = 1.0;
			}
			else
			{
				render->camera.Xiw[i][j] = 0.0;
			}
		}
	}
	/*According to prof, which is SideVec = ViewDir x WorldUp, and UpVector = SideVec x ViewDir*/
	cl = sqrt((render->camera.lookat[X]-render->camera.position[X])*(render->camera.lookat[X]-render->camera.position[X])
					+(render->camera.lookat[Y]-render->camera.position[Y])*(render->camera.lookat[Y]-render->camera.position[Y])
					+(render->camera.lookat[Z]-render->camera.position[Z])*(render->camera.lookat[Z]-render->camera.position[Z]));
	z[X] = (render->camera.lookat[X]-render->camera.position[X])/cl;
	z[Y] = (render->camera.lookat[Y]-render->camera.position[Y])/cl;
	z[Z] = (render->camera.lookat[Z]-render->camera.position[Z])/cl;				/*Init Z*/
	norm = sqrt((z[Z]*render->camera.worldup[Y]-z[Y]*render->camera.worldup[Z])*(z[Z]*render->camera.worldup[Y]-z[Y]*render->camera.worldup[Z])
				+(z[X]*render->camera.worldup[Z]-z[Z]*render->camera.worldup[X])*(z[X]*render->camera.worldup[Z]-z[Z]*render->camera.worldup[X])
				+(z[Y]*render->camera.worldup[X]-z[X]*render->camera.worldup[Y])*(z[Y]*render->camera.worldup[X]-z[X]*render->camera.worldup[Y]));
	x[X] = -(z[Y]*render->camera.worldup[Z]-z[Z]*render->camera.worldup[Y])/norm;
	x[Y] = -(z[Z]*render->camera.worldup[X]-z[X]*render->camera.worldup[Z])/norm;
	x[Z] = -(z[X]*render->camera.worldup[Y]-z[Y]*render->camera.worldup[X])/norm;	/*Init X*/
	norm = sqrt((x[Z]*z[Y]-x[Y]*z[Z])*(x[Z]*z[Y]-x[Y]*z[Z])+(x[X]*z[Z]-x[Z]*z[X])*(x[X]*z[Z]-x[Z]*z[X])+(x[Y]*z[X]-x[X]*z[Y])*(x[Y]*z[X]-x[X]*z[Y]));
	y[X] = -(x[Y]*z[Z]-x[Z]*z[Y])/norm;
	y[Y] = -(x[Z]*z[X]-x[X]*z[Z])/norm;
	y[Z] = -(x[X]*z[Y]-x[Y]*z[X])/norm;												/*Init Y*/
	/*******************************************************************************************/

	/*Compute Xiw*/
	render->camera.Xiw[0][0] = x[X];
	render->camera.Xiw[0][1] = x[Y];
	render->camera.Xiw[0][2] = x[Z];
	render->camera.Xiw[1][0] = y[X];
	render->camera.Xiw[1][1] = y[Y];
	render->camera.Xiw[1][2] = y[Z];
	render->camera.Xiw[2][0] = z[X];
	render->camera.Xiw[2][1] = z[Y];
	render->camera.Xiw[2][2] = z[Z];
	render->camera.Xiw[0][3] = -(x[X]*render->camera.position[X]+x[Y]*render->camera.position[Y]+x[Z]*render->camera.position[Z]);
	render->camera.Xiw[1][3] = -(y[X]*render->camera.position[X]+y[Y]*render->camera.position[Y]+y[Z]*render->camera.position[Z]);
	render->camera.Xiw[2][3] = -(z[X]*render->camera.position[X]+z[Y]*render->camera.position[Y]+z[Z]*render->camera.position[Z]);
	/*********************/

	/*Update Xsp matrix for change of FOV*/
	render->Xsp[2][2] = INT_MAX*tan((render->camera.FOV/2)*PI/180);
	/*************************************/

	/*Initialize Ximage matrix stack*/
	for(i=0;i<MATLEVELS;i++)
	{
		for(j=0;j<4;j++)
		{
			for(k=0;k<4;k++)
			{
				if(k==j)
				{
					render->Ximage[i][j][k] = 1.0;
				}
				else
				{
					render->Ximage[i][j][k] = 0.0;
				}
			}
		}
	}
	/********************************/

	/*Initialize Ximage matrix stack*/
	GzPushMatrix(render,render->Xsp);			/*Push Xsp*/
	GzPushMatrix(render,render->camera.Xpi);	/*Push Xpi*/
	GzPushMatrix(render,render->camera.Xiw);	/*Push Xiw*/
	/********************************/

	render->open = render->display->open;	/*Open render*/
	return GZ_SUCCESS;
}

int GzPutCamera(GzRender *render, GzCamera *camera)
{
/*
- overwrite renderer camera structure with new camera definition
*/
	render->camera.FOV = camera->FOV;					/*Overwrite FOV*/
	render->camera.lookat[X] = camera->lookat[X];		/*Overwrite lookat position*/
	render->camera.lookat[Y] = camera->lookat[Y];
	render->camera.lookat[Z] = camera->lookat[Z];
	render->camera.position[X] = camera->position[X];	/*Overwrite world position*/
	render->camera.position[Y] = camera->position[Y];
	render->camera.position[Z] = camera->position[Z];
	render->camera.worldup[X] = camera->worldup[X];		/*Overwrite world up vector*/
	render->camera.worldup[Y] = camera->worldup[Y];
	render->camera.worldup[Z] = camera->worldup[Z];

	return GZ_SUCCESS;	
}

int GzPushMatrix(GzRender *render, GzMatrix	matrix)
{
/*
- push a matrix onto the Ximage stack
- check for stack overflow
*/
	int i,j;
	GzMatrix dummy=
	{ 
		1.0,	0.0,	0.0,	0.0, 
		0.0,	1.0,	0.0,	0.0, 
		0.0,	0.0,	1.0,	0.0, 
		0.0,	0.0,	0.0,	1.0 
	};

	/*Push matrix to Ximage stack*/
	if(render->matlevel==0)
	{
		for(i=0;i<4;i++)
		{
			for(j=0;j<4;j++)
			{
				render->Ximage[render->matlevel][i][j]=matrix[i][j];
			}
		}
	}
	else
	{
		MatrixMultiply(render->Ximage[render->matlevel-1],matrix,render->Ximage[render->matlevel]);
	}
	/*****************************/

	/*Push matrix to Xnorm stack*/
	if((render->matlevel==0)||(render->matlevel==1))//Push Xsp and Xpi as identity matrix
	{
		for(i=0;i<4;i++)
		{
			for(j=0;j<4;j++)
			{
				render->Xnorm[render->matlevel][i][j]=dummy[i][j];
			}
		}
	}
	else
	{
		CopyMatrix(dummy,matrix);	/*Prevent to modify the original matrix*/
		dummy[0][3]=0.0;			/*Strip translation*/
		dummy[1][3]=0.0;
		dummy[2][3]=0.0;
		Unitary(dummy);				/*Unitary rotation, strip scaling*/
		MatrixMultiply(render->Xnorm[render->matlevel-1],dummy,render->Xnorm[render->matlevel]);
	}
	/****************************/

	render->matlevel++;
	if(render->matlevel>=MATLEVELS)
	{
		return GZ_FAILURE;
	}

	return GZ_SUCCESS;
}

int GzPopMatrix(GzRender *render)
{
/*
- pop a matrix off the Ximage stack
- check for stack underflow
*/
	/*Matlevel minus one first, then pop out the matrix*/
	render->matlevel--;
	if(render->matlevel<0)
	{
		return GZ_FAILURE;
	}

	return GZ_SUCCESS;
}

int GzPutAttribute(GzRender	*render, int numAttributes, GzToken	*nameList, GzPointer *valueList) /* void** valuelist */
{
/*
- set renderer attribute states (e.g.: GZ_RGB_COLOR default color)
- later set shaders, interpolaters, texture maps, and lights
*/
	for(int i=0;i<numAttributes;i++)
	{
		if(nameList[i]==GZ_RGB_COLOR)//Set flat color
		{
			float* color=(float*)valueList[i];
			render->flatcolor[RED]=color[RED];
			render->flatcolor[GREEN]=color[GREEN];
			render->flatcolor[BLUE]=color[BLUE];
		}
		else if(nameList[i]==GZ_INTERPOLATE)//Set interpolater mode
		{
			int* mode=(int*)valueList[i];
			render->interp_mode=*mode;
		}
		else if(nameList[i]==GZ_SHADER)//Set shader mode
		{
		}
		else if(nameList[i]==GZ_DIRECTIONAL_LIGHT)//Set directional light
		{
			GzLight* dLight=(GzLight*)valueList[i];
			render->lights[render->numlights].color[RED]=dLight->color[RED];
			render->lights[render->numlights].color[BLUE]=dLight->color[BLUE];
			render->lights[render->numlights].color[GREEN]=dLight->color[GREEN];
			render->lights[render->numlights].direction[X]=dLight->direction[X];
			render->lights[render->numlights].direction[Y]=dLight->direction[Y];
			render->lights[render->numlights].direction[Z]=dLight->direction[Z];
			render->numlights++;
		}
		else if(nameList[i]==GZ_AMBIENT_LIGHT)//Set ambient light
		{
			GzLight* aLight=(GzLight*)valueList[i];
			render->ambientlight.color[RED]=aLight->color[RED];
			render->ambientlight.color[BLUE]=aLight->color[BLUE];
			render->ambientlight.color[GREEN]=aLight->color[GREEN];
			render->ambientlight.direction[X]=aLight->direction[X];
			render->ambientlight.direction[Y]=aLight->direction[Y];
			render->ambientlight.direction[Z]=aLight->direction[Z];
		}
		else if(nameList[i]==GZ_AMBIENT_COEFFICIENT)//Set ambient coefficient, which is Ka in the render
		{
			float* color=(float*)valueList[i];
			render->Ka[RED]=color[RED];
			render->Ka[GREEN]=color[GREEN];
			render->Ka[BLUE]=color[BLUE];
		}
		else if(nameList[i]==GZ_DIFFUSE_COEFFICIENT)//Set diffuse coefficient, which is Kd in the render
		{
			float* color=(float*)valueList[i];
			render->Kd[RED]=color[RED];
			render->Kd[GREEN]=color[GREEN];
			render->Kd[BLUE]=color[BLUE];
		}
		else if(nameList[i]==GZ_SPECULAR_COEFFICIENT)//Set specular coefficient, which is Ks in the render
		{
			float* color=(float*)valueList[i];
			render->Ks[RED]=color[RED];
			render->Ks[GREEN]=color[GREEN];
			render->Ks[BLUE]=color[BLUE];
		}
		else if(nameList[i]==GZ_DISTRIBUTION_COEFFICIENT)//Set specular power
		{
			float* spec=(float*)valueList[i];
			render->spec=*spec;
		}
		else if(nameList[i]==GZ_AASHIFTX)
		{
			float* shiftx=(float*)valueList[i];
			render->aashiftx=*shiftx;
		}
		else if(nameList[i]==GZ_AASHIFTY)
		{
			float* shifty=(float*)valueList[i];
			render->aashifty=*shifty;
		}
	}
	
	return GZ_SUCCESS;
}

int GzPutTriangle(GzRender	*render, int numParts, GzToken *nameList, GzPointer *valueList)//
/* numParts : how many names and values */
{
/*  
- pass in a triangle description with tokens and values corresponding to 
      GZ_POSITION:3 vert positions in model space /
- Xform positions of verts  /
- Clip - just discard any triangle with verts behind view plane /
       - test for triangles with all three verts off-screen /
- invoke triangle rasterizer  /
*/
	float delta;
	double z;
	int i;
	float tmp[4][1]={1.0,1.0,1.0,1.0};
	float result[4][1]={1.0,1.0,1.0,1.0};
	GzCoord vertex[3];
	GzCoord normal[3];
	GzColor color[3];
	GzIntensity oldr,oldg,oldb,olda;
	GzDepth oldZ;
	//initial edge DDAs and span DDAs
	EdgeDDA* e1=(EdgeDDA*) operator new (sizeof(EdgeDDA));
	EdgeDDA* e2=(EdgeDDA*) operator new (sizeof(EdgeDDA));
	EdgeDDA* e3=(EdgeDDA*) operator new (sizeof(EdgeDDA));
	SpanDDA* s=(SpanDDA*) operator new (sizeof(SpanDDA));
	short red;
	short green;
	short blue;

	//HW2 w/ Xform and Clip
	for(i=0;i<numParts;i++)
	{
		if(nameList[i]==GZ_NULL_TOKEN)
		{
		}
		else if(nameList[i]==GZ_POSITION)//Set position triangle
		{
			float* temp=(float*)valueList[i];
			//get the location of each point
			vertex[0][X]=temp[0];
			vertex[0][Y]=temp[1];
			vertex[0][Z]=temp[2];
			vertex[1][X]=temp[3];
			vertex[1][Y]=temp[4];
			vertex[1][Z]=temp[5];
			vertex[2][X]=temp[6];
			vertex[2][Y]=temp[7];
			vertex[2][Z]=temp[8];
		}
		else if(nameList[i]==GZ_NORMAL)//Set normals of each verts
		{
			float* temp=(float*)valueList[i];
			//get the normals of each point
			normal[0][X]=temp[0];
			normal[0][Y]=temp[1];
			normal[0][Z]=temp[2];
			normal[1][X]=temp[3];
			normal[1][Y]=temp[4];
			normal[1][Z]=temp[5];
			normal[2][X]=temp[6];
			normal[2][Y]=temp[7];
			normal[2][Z]=temp[8];
		}
	}
	/*For flat shading*/
	shade(normal[0],render->flatcolor);
	/******************/

	/*Need to apply matrix here*/
	for(i=0;i<3;i++)
	{
		/*Apply the Ximage stack to vertex*/
		result[0][0]=1.0;
		result[1][0]=1.0;
		result[2][0]=1.0;
		result[3][0]=1.0;
		tmp[0][0]=vertex[i][X];
		tmp[1][0]=vertex[i][Y];
		tmp[2][0]=vertex[i][Z];
		tmp[3][0]=1.0;
		MatrixVectorMultiply(render->Ximage[render->matlevel-1],tmp,result);
		vertex[i][X]=result[0][0]/result[3][0];
		vertex[i][Y]=result[1][0]/result[3][0];
		vertex[i][Z]=result[2][0]/result[3][0];
		/*Apply the Xnorm stack to normal*/
		result[0][0]=1.0;
		result[1][0]=1.0;
		result[2][0]=1.0;
		result[3][0]=1.0;
		tmp[0][0]=normal[i][X];
		tmp[1][0]=normal[i][Y];
		tmp[2][0]=normal[i][Z];
		tmp[3][0]=1.0;
		MatrixVectorMultiply(render->Xnorm[render->matlevel-1],tmp,result);
		normal[i][X]=result[0][0]/result[3][0];
		normal[i][Y]=result[1][0]/result[3][0];
		normal[i][Z]=result[2][0]/result[3][0];
	}
	/***************************/
	
	/*Anti-Aliasing*/
	for(i=0;i<3;i++)
	{
		vertex[i][X]+=render->aashiftx;
		vertex[i][Y]+=render->aashifty;
	}
	/**************/

	/*Discard any triangle with verts behind view plane*/
	if((vertex[0][Z]<0)||(vertex[1][Z]<0)||(vertex[2][Z]<0))
	{
		return GZ_SUCCESS;
	}
	/***************************************************/

	/*Initialize the color*/
	color[0][RED]=0.0;
	color[0][GREEN]=0.0;
	color[0][BLUE]=0.0;
	color[1][RED]=0.0;
	color[1][GREEN]=0.0;
	color[1][BLUE]=0.0;
	color[2][RED]=0.0;
	color[2][GREEN]=0.0;
	color[2][BLUE]=0.0;
	/************************************************************/

	if(render->interp_mode==GZ_COLOR)
	{
		/*Start shading by Gouraud*/
		/*Evaluate the color of each vert by using lighting equation*/
		ComputeColor(render,3,normal,color);
		/************************************************************/

		/*Need to interpolate color for each pixel*/
		/*Rasterization start*/
		YSort(e1,e2,e3,vertex,normal,color);//sort verts by Y
		LRTest(e1,e2,e3);//setup edge DDAs and sort them by L or R
		/**********************************/
		
		if((e1->lr!=T)&&(e2->lr!=B))
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e1);
			AdvanceYCeiling(e3);
			
			Rasterization(render,vertex,normal,color,e1,e3,render->interp_mode);

			AdvanceYCeiling(e2);

			Rasterization(render,vertex,normal,color,e2,e3,render->interp_mode);
		}
		else if(e1->lr==T)
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e2);
			AdvanceYCeiling(e3);

			Rasterization(render,vertex,normal,color,e2,e3,render->interp_mode);
		}
		else if(e2->lr==B)
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e1);
			AdvanceYCeiling(e3);

			Rasterization(render,vertex,normal,color,e1,e3,render->interp_mode);
		}
		/******************************************/
	}
	else if(render->interp_mode==GZ_NORMALS)
	{
		/*Start shading by Phong*/
		/************************/

		/*Need to interpolate normal for each pixel, then compute the color*/
		/*Rasterization start*/
		YSort(e1,e2,e3,vertex,normal,color);//sort verts by Y
		LRTest(e1,e2,e3);//setup edge DDAs and sort them by L or R
		/**********************************/
		
		if((e1->lr!=T)&&(e2->lr!=B))
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e1);
			AdvanceYCeiling(e3);
			
			Rasterization(render,vertex,normal,color,e1,e3,render->interp_mode);

			AdvanceYCeiling(e2);

			Rasterization(render,vertex,normal,color,e2,e3,render->interp_mode);
		}
		else if(e1->lr==T)
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e2);
			AdvanceYCeiling(e3);

			Rasterization(render,vertex,normal,color,e2,e3,render->interp_mode);
		}
		else if(e2->lr==B)
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e1);
			AdvanceYCeiling(e3);

			Rasterization(render,vertex,normal,color,e1,e3,render->interp_mode);
		}
		/******************************************/
	}
	else if(render->interp_mode==GZ_NONE)
	{
		/*Rasterization start*/
		//YSortForGZNone(e1,e2,e3,vertex,normal);//sort verts by Y
		YSort(e1,e2,e3,vertex,normal,color);
		LRTest(e1,e2,e3);//setup edge DDAs and sort them by L or R
		/**********************************/
		
		if((e1->lr!=T)&&(e2->lr!=B))
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e1);
			AdvanceYCeiling(e3);
			
			Rasterization(render,vertex,normal,color,e1,e3,render->interp_mode);

			AdvanceYCeiling(e2);

			Rasterization(render,vertex,normal,color,e2,e3,render->interp_mode);
		}
		else if(e1->lr==T)
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e2);
			AdvanceYCeiling(e3);

			Rasterization(render,vertex,normal,color,e2,e3,render->interp_mode);
		}
		else if(e2->lr==B)
		{
			//Advance and DDA current positions to top y-scan line (ceiling)
			AdvanceYCeiling(e1);
			AdvanceYCeiling(e3);

			Rasterization(render,vertex,normal,color,e1,e3,render->interp_mode);
		}
	}

	/*Rasterization end*/
	delete(e1);
	delete(e2);
	delete(e3);
	delete(s);

	return GZ_SUCCESS;
}

/* NOT part of API - just for general assistance */

short ctoi(float color)		/* convert float color to GzIntensity short */
{
  return(short)((int)(color * ((1 << 12) - 1)));
}

void YSort(EdgeDDA *e1,EdgeDDA *e2,EdgeDDA *e3,GzCoord *vertex,GzCoord *normal,GzColor *color)
{//Sort verts by Y from low to high
	GzCoord temp;
	GzColor ctemp;
	temp[X]=vertex[0][X];
	temp[Y]=vertex[0][Y];
	temp[Z]=vertex[0][Z];
	for(int i=0;i<3;i++)
	{
		for(int j=i;j<3;j++)
		{
			if(vertex[j][Y]<vertex[i][Y])
			{
				/*Change verts*/
				temp[X]=vertex[i][X];
				temp[Y]=vertex[i][Y];
				temp[Z]=vertex[i][Z];
				vertex[i][X]=vertex[j][X];
				vertex[i][Y]=vertex[j][Y];
				vertex[i][Z]=vertex[j][Z];
				vertex[j][X]=temp[X];
				vertex[j][Y]=temp[Y];
				vertex[j][Z]=temp[Z];
				/*Change the noraml of the verts*/
				temp[X]=normal[i][X];
				temp[Y]=normal[i][Y];
				temp[Z]=normal[i][Z];
				normal[i][X]=normal[j][X];
				normal[i][Y]=normal[j][Y];
				normal[i][Z]=normal[j][Z];
				normal[j][X]=temp[X];
				normal[j][Y]=temp[Y];
				normal[j][Z]=temp[Z];
				/*Change the color of the verts*/
				ctemp[RED]=color[i][RED];
				ctemp[GREEN]=color[i][GREEN];
				ctemp[BLUE]=color[i][BLUE];
				color[i][RED]=color[j][RED];
				color[i][GREEN]=color[j][GREEN];
				color[i][BLUE]=color[j][BLUE];
				color[j][RED]=ctemp[RED];
				color[j][GREEN]=ctemp[GREEN];
				color[j][BLUE]=ctemp[BLUE];
			}
			else if(vertex[j][Y]==vertex[i][Y])
			{
				if(vertex[j][X]<=vertex[i][X])
				{
					/*Change verts*/
					temp[X]=vertex[i][X];
					temp[Y]=vertex[i][Y];
					temp[Z]=vertex[i][Z];
					vertex[i][X]=vertex[j][X];
					vertex[i][Y]=vertex[j][Y];
					vertex[i][Z]=vertex[j][Z];
					vertex[j][X]=temp[X];
					vertex[j][Y]=temp[Y];
					vertex[j][Z]=temp[Z];
					/*Change the noraml of the verts*/
					temp[X]=normal[i][X];
					temp[Y]=normal[i][Y];
					temp[Z]=normal[i][Z];
					normal[i][X]=normal[j][X];
					normal[i][Y]=normal[j][Y];
					normal[i][Z]=normal[j][Z];
					normal[j][X]=temp[X];
					normal[j][Y]=temp[Y];
					normal[j][Z]=temp[Z];
					/*Change the color of the verts*/
					ctemp[RED]=color[i][RED];
					ctemp[GREEN]=color[i][GREEN];
					ctemp[BLUE]=color[i][BLUE];
					color[i][RED]=color[j][RED];
					color[i][GREEN]=color[j][GREEN];
					color[i][BLUE]=color[j][BLUE];
					color[j][RED]=ctemp[RED];
					color[j][GREEN]=ctemp[GREEN];
					color[j][BLUE]=ctemp[BLUE];
				}
			}
		}
	}
	//Initial edge DDAs
	if(vertex[0][Y]==vertex[1][Y])
	{//if the top edge is horizontal
		InitEdgeDDAT(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);

	}
	else if(vertex[1][Y]==vertex[2][Y])
	{//if the bottom edge is horizontal
		InitEdgeDDAB(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);
	}
	else
	{
		InitEdgeDDA(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);
	}
}

void YSortForGZNone(EdgeDDA *e1,EdgeDDA *e2,EdgeDDA *e3,GzCoord *vertex,GzCoord *normal)
{//Sort verts by Y from low to high
	GzCoord temp;
	temp[X]=vertex[0][X];
	temp[Y]=vertex[0][Y];
	temp[Z]=vertex[0][Z];
	for(int i=0;i<3;i++)
	{
		for(int j=i;j<3;j++)
		{
			if(vertex[j][Y]<vertex[i][Y])
			{
				/*Change verts*/
				temp[X]=vertex[i][X];
				temp[Y]=vertex[i][Y];
				temp[Z]=vertex[i][Z];
				vertex[i][X]=vertex[j][X];
				vertex[i][Y]=vertex[j][Y];
				vertex[i][Z]=vertex[j][Z];
				vertex[j][X]=temp[X];
				vertex[j][Y]=temp[Y];
				vertex[j][Z]=temp[Z];
				/*Change the noraml of the verts*/
				temp[X]=normal[i][X];
				temp[Y]=normal[i][Y];
				temp[Z]=normal[i][Z];
				normal[i][X]=normal[j][X];
				normal[i][Y]=normal[j][Y];
				normal[i][Z]=normal[j][Z];
				normal[j][X]=temp[X];
				normal[j][Y]=temp[Y];
				normal[j][Z]=temp[Z];
			}
			else if(vertex[j][Y]==vertex[i][Y])
			{
				if(vertex[j][X]<=vertex[i][X])
				{
					/*Change verts*/
					temp[X]=vertex[i][X];
					temp[Y]=vertex[i][Y];
					temp[Z]=vertex[i][Z];
					vertex[i][X]=vertex[j][X];
					vertex[i][Y]=vertex[j][Y];
					vertex[i][Z]=vertex[j][Z];
					vertex[j][X]=temp[X];
					vertex[j][Y]=temp[Y];
					vertex[j][Z]=temp[Z];
					/*Change the noraml of the verts*/
					temp[X]=normal[i][X];
					temp[Y]=normal[i][Y];
					temp[Z]=normal[i][Z];
					normal[i][X]=normal[j][X];
					normal[i][Y]=normal[j][Y];
					normal[i][Z]=normal[j][Z];
					normal[j][X]=temp[X];
					normal[j][Y]=temp[Y];
					normal[j][Z]=temp[Z];
				}
			}
		}
	}
	//Initial edge DDAs
	if(vertex[0][Y]==vertex[1][Y])
	{//if the top edge is horizontal
		InitEdgeDDAT(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);

	}
	else if(vertex[1][Y]==vertex[2][Y])
	{//if the bottom edge is horizontal
		InitEdgeDDAB(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);
	}
	else
	{
		InitEdgeDDA(e1,vertex[X],vertex[Y]);
		InitEdgeDDA(e2,vertex[Y],vertex[Z]);
		InitEdgeDDA(e3,vertex[X],vertex[Z]);
	}
}

void InitEdgeDDA(EdgeDDA *e,float *v1,float *v2)
{//Initial edge DDA
	e->vStart[X]=v1[X];
	e->vStart[Y]=v1[Y];
	e->vStart[Z]=v1[Z];
	e->vEnd[X]=v2[X];
	e->vEnd[Y]=v2[Y];
	e->vEnd[Z]=v2[Z];
	e->vCurrent[X]=v1[X];
	e->vCurrent[Y]=v1[Y];
	e->vCurrent[Z]=v1[Z];
	e->slopex=(v2[X]-v1[X])/(v2[Y]-v1[Y]);
	e->slopez=(v2[Z]-v1[Z])/(v2[Y]-v1[Y]);
	e->lr=L;
}

void InitEdgeDDAT(EdgeDDA *e,float *v1,float *v2)
{//Initial horizontal top edge DDA 
	e->vStart[X]=v1[X];
	e->vStart[Y]=v1[Y];
	e->vStart[Z]=v1[Z];
	e->vEnd[X]=v2[X];
	e->vEnd[Y]=v2[Y];
	e->vEnd[Z]=v2[Z];
	e->vCurrent[X]=v1[X];
	e->vCurrent[Y]=v1[Y];
	e->vCurrent[Z]=v1[Z];
	e->slopex=0;
	e->slopez=0;
	e->lr=T;
}

void InitEdgeDDAB(EdgeDDA *e,float *v1,float *v2)
{//Initial horizontal bottom edge DDA
	e->vStart[X]=v1[X];
	e->vStart[Y]=v1[Y];
	e->vStart[Z]=v1[Z];
	e->vEnd[X]=v2[X];
	e->vEnd[Y]=v2[Y];
	e->vEnd[Z]=v2[Z];
	e->vCurrent[X]=v1[X];
	e->vCurrent[Y]=v1[Y];
	e->vCurrent[Z]=v1[Z];
	e->slopex=0;
	e->slopez=0;
	e->lr=B;
}

void LRTest(EdgeDDA *e1,EdgeDDA *e2,EdgeDDA *e3)
{//Sort edge DDA by L or R
	if(e1->lr==T)
	{//if top edge is horizontal
		e3->lr=L;
		e2->lr=R;
	}
	else if(e2->lr==B)
	{//if bottom edge is horizontal
		e1->lr=L;
		e3->lr=R;
	}
	else
	{
		if(e1->slopex<e3->slopex)
		{
			e3->lr=R;
		}
		else
		{
			e1->lr=R;
			e2->lr=R;
		}//sort L/R by slope
	}
}

void InitSpanDDA(SpanDDA *s,EdgeDDA *e1,EdgeDDA *e2)
{//Initial span DDA
	if(e1->lr==L)
	{//if edge1 is the left edge
		if(e1->vCurrent[X]<0)
		{//if the left end is out of the bound, start the scan line from x=0
			s->xStart=0;
			s->slopez=(e2->vCurrent[Z]-e1->vCurrent[Z])/(e2->vCurrent[X]-e1->vCurrent[X]);
			s->zStart=e1->vCurrent[Z]+s->slopez*(0-e1->vCurrent[X]);
			s->xEnd=e2->vCurrent[X];
			s->zEnd=e2->vCurrent[Z];
			s->xCurrent=s->xStart;
			s->zCurrent=s->zStart;
		}
		else
		{
			s->xStart=e1->vCurrent[X];
			s->zStart=e1->vCurrent[Z];
			s->xEnd=e2->vCurrent[X];
			s->zEnd=e2->vCurrent[Z];
			s->xCurrent=e1->vCurrent[X];
			s->zCurrent=e1->vCurrent[Z];
			s->slopez=(s->zEnd-s->zStart)/(s->xEnd-s->xStart);
		}
	}
	else if(e2->lr==L)
	{//if edge2 is left edge
		if(e2->vCurrent[X]<0)
		{//if the left end is out of the bound, start the scan line from x=0
			s->xStart=0;
			s->slopez=(e1->vCurrent[Z]-e2->vCurrent[Z])/(e1->vCurrent[X]-e2->vCurrent[X]);
			s->zStart=e2->vCurrent[Z]+s->slopez*(0-e2->vCurrent[X]);
			s->xEnd=e1->vCurrent[X];
			s->zEnd=e1->vCurrent[Z];
			s->xCurrent=s->xStart;
			s->zCurrent=s->zStart;
		}
		else
		{
			s->xStart=e2->vCurrent[X];
			s->zStart=e2->vCurrent[Z];
			s->xEnd=e1->vCurrent[X];
			s->zEnd=e1->vCurrent[Z];
			s->xCurrent=e2->vCurrent[X];
			s->zCurrent=e2->vCurrent[Z];
			s->slopez=(s->zEnd-s->zStart)/(s->xEnd-s->xStart);
		}
	}
}

void AdvanceYCeiling(EdgeDDA *e)
{//Advance and DDA current positions to top y-scan line (ceiling)
	float delta;
	delta=(int)e->vStart[Y]+1-e->vStart[Y];
	if(delta!=1)
	{
		e->vCurrent[X]=e->vCurrent[X]+e->slopex*delta;
		e->vCurrent[Y]=e->vCurrent[Y]+delta;
		e->vCurrent[Z]=e->vCurrent[Z]+e->slopez*delta;
	}
}

void MatrixMultiply (float (*a)[4], float (*b)[4], float (*c)[4])
{//Multiplication between two 4x4 matrixes, matrix a and b are the source and matrix c is the target
	int i,j,k;
	float tmp;

	for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			tmp=0.0;
			for(k=0;k<4;k++)
			{
				tmp+=a[i][k]*b[k][j];
			}
			c[i][j]=tmp;
		}
	}
	return;
}

void MatrixVectorMultiply (float (*a)[4], float (*b)[1], float (*c)[1])
{//Multiplication between one 4x4 matrix and one 4x1 matrix, and matrix a is 4x4, matrix b is 4x1, and the target matrix c is 4x1
	int i,j,k;
	float tmp;

	for(i=0;i<4;i++)
	{
		for(j=0;j<1;j++)
		{
			tmp=0.0;
			for(k=0;k<4;k++)
			{
				tmp+=a[i][k]*b[k][j];
			}
			c[i][j]=tmp;
		}
	}
	return;
}

void CopyMatrix(float (*target)[4],float (*source)[4])
{//Copy one 4x4 matrix into another 4x4 matrix
	int i,j;

	for(i=0;i<4;i++)
	{
		for(j=0;j<4;j++)
		{
			target[i][j]=source[i][j];
		}
	}
}

//HW4
int Unitary(GzMatrix matrix)
{//Unitary rotation for each matrix
	float k=0;
	int i,j;
	/*Calculate the scale factor*/
	for(i=0;i<3;i++)
	{
		k+=matrix[i][0]*matrix[i][0]+matrix[i][1]*matrix[i][1]+matrix[i][2]*matrix[i][2];
	}
	k=sqrt(k);
	/****************************/
	/*Apply the factor to all elements*/
	for(i=0;i<3;i++)
	{
		for(j=0;j<3;j++)
		{
			matrix[i][j]=matrix[i][j]/k;
		}
	}
	/**********************************/
	return GZ_SUCCESS;
}

void GouraudInterpolation(float (*color)[3],float (*vertex)[3],short &red,short &green,short &blue,float x,float y)
{//Gouraud interpolation
	float lambda1,lambda2,lambda3;
	float r,g,b;

	/*Using Barycentric for interpolation*/
	lambda1=((vertex[1][Y]-vertex[2][Y])*(x-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(y-vertex[2][Y]))/((vertex[1][Y]-vertex[2][Y])*(vertex[0][X]-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(vertex[0][Y]-vertex[2][Y]));
	lambda2=((vertex[2][Y]-vertex[0][Y])*(x-vertex[2][X])+(vertex[0][X]-vertex[2][X])*(y-vertex[2][Y]))/((vertex[1][Y]-vertex[2][Y])*(vertex[0][X]-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(vertex[0][Y]-vertex[2][Y]));
	lambda3=1-lambda1-lambda2;
	/************************************/
	/*Interpolate color for every pixel*/
	r=lambda1*color[0][RED]+lambda2*color[1][RED]+lambda3*color[2][RED];
	g=lambda1*color[0][GREEN]+lambda2*color[1][GREEN]+lambda3*color[2][GREEN];
	b=lambda1*color[0][BLUE]+lambda2*color[1][BLUE]+lambda3*color[2][BLUE];
	/***********************************/
	/*Check for color overflow*/
	if(r>1)
	{
		r=1;
	}
	if(g>1)
	{
		g=1;
	}
	if(b>1)
	{
		b=1;
	}
	/**************************/
	/*Convert float value to 12-bit*/
	red=ctoi(r);
	green=ctoi(g);
	blue=ctoi(b);
	/*******************************/
}

void PhongInterpolation(GzRender *render,GzCoord *vertex,GzCoord *normal,short &red,short &green,short &blue,float x,float y)
{//Phong Interpolation for pixel color
	float lambda1,lambda2,lambda3;
	float r,g,b;
	GzColor color={0.0,0.0,0.0};
	GzCoord normalCurrent;

	/*Using Barycentric for interpolation*/
	lambda1=((vertex[1][Y]-vertex[2][Y])*(x-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(y-vertex[2][Y]))/((vertex[1][Y]-vertex[2][Y])*(vertex[0][X]-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(vertex[0][Y]-vertex[2][Y]));
	lambda2=((vertex[2][Y]-vertex[0][Y])*(x-vertex[2][X])+(vertex[0][X]-vertex[2][X])*(y-vertex[2][Y]))/((vertex[1][Y]-vertex[2][Y])*(vertex[0][X]-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(vertex[0][Y]-vertex[2][Y]));
	lambda3=1-lambda1-lambda2;
	/************************************/
	/*Interpolate normal for every pixel*/
	normalCurrent[X]=lambda1*normal[0][X]+lambda2*normal[1][X]+lambda3*normal[2][X];
	normalCurrent[Y]=lambda1*normal[0][Y]+lambda2*normal[1][Y]+lambda3*normal[2][Y];
	normalCurrent[Z]=lambda1*normal[0][Z]+lambda2*normal[1][Z]+lambda3*normal[2][Z];
	/************************************/

	ComputeColor(render,1,&normalCurrent,&color);//Compute color according to the own normal

	r=color[RED];
	g=color[GREEN];
	b=color[BLUE];
	/*Check for color overflow*/
	if(r>1)
	{
		r=1;
	}
	if(g>1)
	{
		g=1;
	}
	if(b>1)
	{
		b=1;
	}
	/**************************/
	/*Convert float value to 12-bit*/
	red=ctoi(r);
	green=ctoi(g);
	blue=ctoi(b);
	/*******************************/
}

double ZInterpolation(GzCoord *vertex,float x,float y)
{//Interpolate Z for every pixel
	double result;
	float lambda1,lambda2,lambda3;

	/*Using Barycentric for interpolation*/
	lambda1=((vertex[1][Y]-vertex[2][Y])*(x-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(y-vertex[2][Y]))/((vertex[1][Y]-vertex[2][Y])*(vertex[0][X]-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(vertex[0][Y]-vertex[2][Y]));
	lambda2=((vertex[2][Y]-vertex[0][Y])*(x-vertex[2][X])+(vertex[0][X]-vertex[2][X])*(y-vertex[2][Y]))/((vertex[1][Y]-vertex[2][Y])*(vertex[0][X]-vertex[2][X])+(vertex[2][X]-vertex[1][X])*(vertex[0][Y]-vertex[2][Y]));
	lambda3=1-lambda1-lambda2;
	/************************************/

	result=lambda1*vertex[0][Z]+lambda2*vertex[1][Z]+lambda3*vertex[2][Z];

	return result;
}

void ComputeColor(GzRender *render,int numverts,GzCoord *normal,GzColor *color)
{//Compute the color by lighting equation
	int i,j;
	GzCoord l,n,r,e;
	float nl,ne,re,sum;

	/*Initialize the vector E by default(0,0,-1)*/
	e[X]=0.0;
	e[Y]=0.0;
	e[Z]=-1.0;
	/********************************************/

	for(i=0;i<numverts;i++)
	{//Calculate for every vertice
		for(j=0;j<render->numlights;j++)
		{//Calculate for each light
			l[X]=render->lights[j].direction[X];
			l[Y]=render->lights[j].direction[Y];
			l[Z]=render->lights[j].direction[Z];
			n[X]=normal[i][X];
			n[Y]=normal[i][Y];
			n[Z]=normal[i][Z];
			/*Normalize the normal of the vertice*/
			sum=sqrt(n[X]*n[X]+n[Y]*n[Y]+n[Z]*n[Z]);
			n[X]=n[X]/sum;
			n[Y]=n[Y]/sum;
			n[Z]=n[Z]/sum;
			nl=DotProduct(n,l);
			/*Calculate vector R*/
			r[X]=2*n[X]*(nl)-l[X];
			r[Y]=2*n[Y]*(nl)-l[Y];
			r[Z]=2*n[Z]*(nl)-l[Z];
			re=DotProduct(r,e);
			ne=DotProduct(n,e);
			if((nl>=0)&&(ne>=0))
			{//Both positive
				if(re<0)
				{//Clamp R dot E to zero
					re=0;
				}
				/*C=Ks*le*(R dot E)^spec+Kd*le*(N dot L) for each directional light*/
				color[i][RED]+=(render->Ks[RED]*render->lights[j].color[RED]*pow(re,render->spec))+(render->Kd[RED]*render->lights[j].color[RED]*nl);
				color[i][GREEN]+=(render->Ks[GREEN]*render->lights[j].color[GREEN]*pow(re,render->spec))+(render->Kd[GREEN]*render->lights[j].color[GREEN]*nl);
				color[i][BLUE]+=(render->Ks[BLUE]*render->lights[j].color[BLUE]*pow(re,render->spec))+(render->Kd[BLUE]*render->lights[j].color[BLUE]*nl);
			}
			else if((nl<0)&&(ne<0))
			{//Both negative, filp the normal
				n[X]=-n[X];
				n[Y]=-n[Y];
				n[Z]=-n[Z];
				nl=DotProduct(n,l);
				r[X]=2*n[X]*(nl)-l[X];
				r[Y]=2*n[Y]*(nl)-l[Y];
				r[Z]=2*n[Z]*(nl)-l[Z];
				re=DotProduct(r,e);
				ne=DotProduct(n,e);
				if(re<0)
				{
					re=0;
				}
				color[i][RED]+=(render->Ks[RED]*render->lights[j].color[RED]*pow(re,render->spec))+(render->Kd[RED]*render->lights[j].color[RED]*nl);
				color[i][GREEN]+=(render->Ks[GREEN]*render->lights[j].color[GREEN]*pow(re,render->spec))+(render->Kd[GREEN]*render->lights[j].color[GREEN]*nl);
				color[i][BLUE]+=(render->Ks[BLUE]*render->lights[j].color[BLUE]*pow(re,render->spec))+(render->Kd[BLUE]*render->lights[j].color[BLUE]*nl);
			}
			else if(((nl>=0)&&(ne<0))||((nl<0)&&(ne>=0)))
			{//Both different, skip
				//return GZ_SUCCESS;
				continue;
			}
		}
		/*C=C+Kd*la*/
		color[i][RED]+=render->Ka[RED]*render->ambientlight.color[RED];
		color[i][GREEN]+=render->Ka[GREEN]*render->ambientlight.color[GREEN];
		color[i][BLUE]+=render->Ka[BLUE]*render->ambientlight.color[BLUE];
	}
}

float DotProduct(GzCoord a,GzCoord b)
{//Calculate the dot product of two vectors
	float result;
	result=a[X]*b[X]+a[Y]*b[Y]+a[Z]*b[Z];
	return result;
}

float VectorLength(GzCoord a)
{//Calculate the length of the vector
	float result;
	result=sqrt(a[X]*a[X]+a[Y]*a[Y]+a[Z]*a[Z]);
	return result;
}

void shade(GzCoord norm, GzCoord color)
{//Calculate the color for flat shading
	GzCoord light;
	float coef;

	light[0] = 0.707f;
	light[1] = 0.5f;
	light[2] = 0.5f;

	coef = light[0]*norm[0] + light[1]*norm[1] + light[2]*norm[2];
	if (coef < 0) 	coef *= -1;

	if (coef > 1.0)	coef = 1.0;
	color[0] = coef*0.95f;
	color[1] = coef*0.65f;
	color[2] = coef*0.88f;
}

void Rasterization(GzRender *render,GzCoord *vertex,GzCoord *normal,GzColor *color,EdgeDDA *e1,EdgeDDA *e2,int interpMode)
{//Rasterize the whole image by using scan-line algorithm
	float delta;
	GzIntensity oldr,oldg,oldb,olda;
	GzDepth oldZ;
	short red,green,blue;
	double z;
	SpanDDA* s=(SpanDDA*) operator new (sizeof(SpanDDA));

	while(e1->vCurrent[Y]<=e1->vEnd[Y])
	{//Start from the lower edge
		InitSpanDDA(s,e1,e2);//Initial span DDA				
		delta=(int)s->xStart+1-s->xStart;
		if(delta!=1)
		{
			s->xCurrent=s->xStart+delta;
			s->zCurrent=s->zStart+s->slopez*delta;
		}
		else
		{
			s->xCurrent=s->xStart;
			s->zCurrent=s->zStart;
		}//Setup span edges and set positions

		while(s->xCurrent<=s->xEnd)
		{
			z=ZInterpolation(vertex,s->xCurrent,e1->vCurrent[Y]);//Interpolate z
			s->zCurrent=z;
			GzGetDisplay(render->display,s->xCurrent,e1->vCurrent[Y],&oldr,&oldg,&oldb,&olda,&oldZ);//Get the z-buffer
			
			if(z<=oldZ)//Check the z with the z-buffer
			{
				if(interpMode==GZ_COLOR)
				{//Gouraud Interpolation
					GouraudInterpolation(color,vertex,red,green,blue,s->xCurrent,e1->vCurrent[Y]);
				}
				else if(interpMode==GZ_NORMALS)
				{//Phong Interpolation
					PhongInterpolation(render,vertex,normal,red,green,blue,s->xCurrent,e1->vCurrent[Y]);
				}
				else if(interpMode==GZ_NONE)
				{//Flat shading
					red=ctoi(render->flatcolor[RED]);
					green=ctoi(render->flatcolor[GREEN]);
					blue=ctoi(render->flatcolor[BLUE]);
				}
				GzPutDisplay(render->display,(int)s->xCurrent,(int)e1->vCurrent[Y],red,green,blue,1,z);//Test interpolate-Z and write color value into FB
			}
			s->xCurrent+=1;//Go to next pixel
			s->zCurrent+=s->slopez;
		}
		//Go to next line
		e1->vCurrent[Y]+=1;
		e1->vCurrent[X]=e1->vCurrent[X]+e1->slopex;
		e1->vCurrent[Z]=e1->vCurrent[Z]+e1->slopez;
		e2->vCurrent[Y]+=1;
		e2->vCurrent[X]=e2->vCurrent[X]+e2->slopex;
		e2->vCurrent[Z]=e2->vCurrent[Z]+e2->slopez;
	}
	delete(s);
	return;
}

int GzAASampleCombination(GzDisplay **display,GzDisplay *result,float (*AAFilter)[3])
{
	int i,j;
	GzIntensity red=0;
	GzIntensity blue=0;
	GzIntensity green=0;

	for(i=0;i<result->xres*result->yres;i++)
	{
		for(j=0;j<AAKERNEL_SIZE;j++)
		{
			red+=display[j]->fbuf[i].red*AAFilter[j][2];
			green+=display[j]->fbuf[i].green*AAFilter[j][2];
			blue+=display[j]->fbuf[i].blue*AAFilter[j][2];
		}
		result->fbuf[i].red=red;
		result->fbuf[i].green=green;
		result->fbuf[i].blue=blue;
		red=0;
		green=0;
		blue=0;
	}
	return GZ_SUCCESS;
}